Represents districts and sites that are listed on the National Register of Historic Places (Federal designation) and the Virginia Landmarks Register (State designation).
For our Fan District analysis we will be working with Civic Associations to get the formal boundary of the Fan District Association.
We then use that boundary to determine Addresses and Parcels in the Fan District Association
Fan District Association
FDA - Key measures
Area
Fan District Association: Area: 24,564,893 sq.ft., 563.93 acres, 0.881 sq.miles
Fan District Association: Area: 2,282,153 sq.m, 228.215 hectares, 2.282 sq.km
Examples of records sharing the same AddressLabel
Total duplicated AddressLabel records: 2
Number of unique duplicated labels: 1
Examples of duplicated AddressLabels:
AddressLabel UnitType UnitValue
100329 801 W Franklin St Room 503 Room 503
100338 801 W Franklin St Room 503 Room 503
Examples of non-unit addresses with unit parts in AddressLabel
🔍 Found 716 rows where AddressLabel ends with a recognizable UnitType + value and UnitType is missing.
🧪 Sample rows:
AddressLabel UnitType UnitValue
63 517 N Arthur Ashe Blvd Unit G9 None None
70 517 N Arthur Ashe Blvd Unit G1 None None
138 517 N Arthur Ashe Blvd Unit G27 None None
149 517 N Arthur Ashe Blvd Unit G15 None None
154 517 N Arthur Ashe Blvd Unit G14 None None
Fix these issues found above.
✅ Extracted UnitType and UnitValue for 716 rows.
Examples of non-unit addreses sharing the same long/lat
🔍 First record from each duplicate (UnitType is None, by lat/lon):
1403 Floyd Ave | 34 records | (37.54774625, -77.45877807)
708 1/2 W Grace St | 2 records | (37.54886389, -77.44932289)
800 W Grace St | 3 records | (37.5494549, -77.45003694)
2601 1/2 W Main St | 2 records | (37.55178826, -77.47347259)
1625 1/2 W Grace St | 3 records | (37.55373633, -77.45790774)
508 1 Allison St | 2 records | (37.55572891, -77.46464775)
2505 C Stuart Ave | 2 records | (37.55583431, -77.46986291)
500 1/2 N Stafford Ave | 2 records | (37.55655244, -77.46784226)
502 1/2 N Stafford Ave | 2 records | (37.55661253, -77.46778968)
2500 1/2 Kensington Ave | 2 records | (37.55715189, -77.46891095)
2502 1/2 Kensington Ave | 2 records | (37.55721386, -77.46903156)
2226 1/2 W Grace St | 2 records | (37.55783315, -77.46423125)
2620 W Grace St | 2 records | (37.56040195, -77.46847792)
2624 W Grace St | 2 records | (37.56043825, -77.46854847)
2628 W Grace St | 2 records | (37.56047908, -77.46862291)
2632 W Grace St | 2 records | (37.56052296, -77.46869797)
Fix: jitter addresses
✅ Jittered coordinates applied to duplicate lat/lon groups.
Example: Unit addresses without non-unit base address
Total unit addresses: 7326
Unit addresses with missing base address: 1
Examples of missing base addresses:
AddressLabel BaseAddressLabel UnitType UnitValue
90017 1501 Hanover Unit A 1501 Hanover Unit A
List of non-unit addresses where baseaddresslabel <> addressLabel
Total base addresses (UnitType is null): 3987
Matching AddressLabel == BaseAddressLabel: 3878
Non-matching AddressLabel != BaseAddressLabel: 109
Examples of mismatched base addresses:
AddressLabel BaseAddressLabel
35 2032 Park Ave Rear 2032 Park Ave
43 1811 Monument Ave Rear 1811 Monument Ave
54 1815 Monument Ave Rear 1815 Monument Ave
64 211 N Arthur Ashe Blvd Rear 211 N Arthur Ashe Blvd
97 1711 Hanover Ave Rear 1711 Hanover Ave
Addresses with the top 20 unit counts.
AddressLabel BaseAddressLabel UnitCount UnitType
77001 710 W Franklin St 710 W Franklin St 981 None
57411 1025 W Grace St 1025 W Grace St 255 None
70222 406 Shafer St 406 Shafer St 155 None
17023 1363 W Broad St 1363 W Broad St 114 None
68256 2501 W Broad St 2501 W Broad St 109 None
13924 801 W Franklin St 801 W Franklin St 94 None
23850 933 W Broad St 933 W Broad St 90 None
69904 1600 Monument Ave 1600 Monument Ave 64 None
43551 900 W Franklin St 900 W Franklin St 63 None
74594 1333 W Broad St RSPS 1333 W Broad St 58 None
23632 1333 W Broad St 1333 W Broad St 58 None
70548 612 W Franklin St 612 W Franklin St 56 None
83188 501 N Allen Ave 501 N Allen Ave 55 None
70067 1030 W Franklin St 1030 W Franklin St 54 None
77910 413 Stuart Cir 413 Stuart Cir 37 None
69576 1005 Grove Ave 1005 Grove Ave 36 None
84307 2621 Stuart Ave 2621 Stuart Ave 36 None
13398 1301 W Broad St 1301 W Broad St 36 None
63125 1630 Monument Ave 1630 Monument Ave 36 None
69481 1108 W Franklin St 1108 W Franklin St 36 None
FDA - Static Map
Below is a static map of the Fan District Association, using data obtained from the City of Richmond GeoHub. The image shows the boundary of the Fan District Association in black. Addresses are identifed by red dots. Parcels are identified as gray regions.
Below is an interactive map of the Fan District Association. You can zoom and scroll.
Make this Notebook Trusted to load map: File -> Trust Notebook
Interactive Map: Heatmap of Addresses and Unitcounts
Make this Notebook Trusted to load map: File -> Trust Notebook
Interactive Map: Parcels and UnitCounts
Make this Notebook Trusted to load map: File -> Trust Notebook
Available columns
Civic_Associations
Property
ID
Name
AdoptionDate
ChangeDate
GlobalID
Shape__Area
Shape__Length
geometry
Addresses_in_fan
Property
AddressId
AddressLabel
BuildingNumber
StreetDirection
StreetName
StreetType
ExtensionWithUnit
UnitType
UnitValue
ZipCode
Mailable
StatePlaneX
StatePlaneY
Latitude
Longitude
geometry
ID
Name
GlobalID
Shape__Area
Shape__Length
BaseAddressLabel
UnitCount
Parcels_in_fan
Property
ParcelID
PIN
CountOfPIN
OwnerName
AsrLocationBldgNo
MailAddress
MailCity
MailState
MailZip
AssessmentDate
LandValue
DwellingValue
TotalValue
LandSqFt
ProvalAsmtNhood
TaxExemptCode
PropertyClassID
PropertyClass
LandUse
Mailable
MaskedOwner
GlobalID_left
geometry
ID
Name
GlobalID_right
Shape__Area
Shape__Length
Source Code
---title: Fan District Association---The City of Richmond maintains an arcgis geodata repository called [Richmond GeoHub](https://richmond-geo-hub-cor.hub.arcgis.com).Data on the geohub are organized into key areas. For our analysis, we'll be using the following data sources.[Addresses](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/674d645c444f4191998f0ebb96e56047_0/explore?location=37.527383%2C-77.493413%2C10.99): All of the official, mapped inventory of all unit and non-unit-based addresses in the City. Includes only active addresses.[Parcels](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/fbfce2aab2a44c05bc0abc2d6ea7e54a_0/explore?location=37.525465%2C-77.493422%2C10.60): City of Richmond property ownership information, mapped by land ownership (parcels).[Civic Associations](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/be39ce592f3e4419babe11d1b967e2f3_0/explore?location=37.528836%2C-77.494197%2C10.96): Represents civic organization boundaries in the city of Richmond, Virginia.[National Historic Districts](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/38bd0df47c6440528c2ef22daaf81883_0/explore?location=37.550339%2C-77.468606%2C14.93): Represents districts and sites that are listed on the National Register of Historic Places (Federal designation) and the Virginia Landmarks Register (State designation).[Neighborhoods](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/7a0ffef23d16461e9728c065f27b2790_0/explore?location=37.525021%2C-77.493427%2C10.73): City of Richmond Neighborhoods.For our Fan District analysis we will be working with [Civic Associations](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/be39ce592f3e4419babe11d1b967e2f3_0/explore?location=37.528836%2C-77.494197%2C10.96)to get the formal boundary of the *Fan District Association*.We then use that boundary to determine [Addresses](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/674d645c444f4191998f0ebb96e56047_0/explore?location=37.527383%2C-77.493413%2C10.99)and [Parcels](https://richmond-geo-hub-cor.hub.arcgis.com/datasets/fbfce2aab2a44c05bc0abc2d6ea7e54a_0/explore?location=37.525465%2C-77.493422%2C10.60)in the *Fan District Association*```{python}import osimport reimport sysimport sqlite3import requestsimport numpy as npimport pandas as pdimport geopandas as gpdimport matplotlib.pyplot as pltfrom matplotlib.patches import Patchfrom matplotlib.lines import Line2Dfrom math import floorfrom loguru import loggersys.path.append("..")from fandu.geo_utils import get_newest_feature_filepd.set_option("display.max_rows", None)geojson_folder ="../precious/"features = ["Addresses","Parcels"]selector ="Civic_Associations"selector_key ="Fan District Association"```# Fan District Association```{python}# Load the neighborhoods GeoJSONdata = {}for feature in [selector] + features: geofile = get_newest_feature_file( geojson_folder, feature )#logger.debug(geofile) data[feature] = gpd.read_file( geofile )for feature in features: data[feature] = data[feature].to_crs( data[selector].crs )# columns to drop:shared_drops = ['CreatedBy','CreatedDate','EditBy','EditDate']drop_columns = {"Civic_Associations" : ['OBJECTID'] + shared_drops,"Addresses" : ['OBJECTID'] + shared_drops,"Parcels" : ['OBJECTID'],}for feature in [selector] + features: data[feature] = data[feature].drop(columns=drop_columns[feature])``````{python}for feature in features: data[selector_key] = data[selector][ data[selector]["Name"] == selector_key ]for feature in features: predicate ="overlaps"if feature=="Neighborhoods"else"within" data[feature+"_in_fan"] = gpd.sjoin(data[feature], data[selector_key], predicate=predicate, how="inner")for feature in features: feature_name = feature+"_in_fan" data[feature_name] = data[feature_name].drop(columns=["index_right","AdoptionDate","ChangeDate"])for feature in features: feature_name = feature+"_in_fan" df = data[feature_name].drop(columns="geometry") df.to_csv(f"{feature_name}.csv", index=False)``````{python}if1:# Connect to or create a database file conn = sqlite3.connect("fda-data.sqlite")for feature in [selector]+features:# Write the DataFrame to a new table (overwrite if it exists) temp_copy = data[feature].copy() temp_copy = temp_copy.drop(columns=['geometry']) temp_copy.to_sql(feature.lower(), conn, if_exists="replace", index=False)# write out _in_fan featuresifnot feature==selector: feature_name = feature +"_in_fan" temp_copy = data[feature_name].copy() temp_copy = temp_copy.drop(columns=['geometry']) temp_copy.to_sql(feature_name.lower(), conn, if_exists="replace", index=False) conn.close()```## FDA - Key measures### Area ```{python}# Get the Shape_Area value for the row where Name == selector_keyprojected = data[selector].to_crs(epsg=2283)# Get Shape_Area for the selected featurearea_sqft = projected.loc[projected["Name"] == selector_key, "geometry"].area.iloc[0]# Convert to acresarea_acres = area_sqft /43560area_sqmi = area_sqft /27_878_400area_sqm = area_sqft *0.09290304area_hectare = area_sqm /10_000area_sqkm = area_sqm /1_000_000print(f"{selector_key}: Area: {area_sqft:,.0f} sq.ft., {area_acres:,.2f} acres, {area_sqmi:,.3f} sq.miles")print(f"{selector_key}: Area: {area_sqm:,.0f} sq.m, {area_hectare:,.3f} hectares, {area_sqkm:,.3f} sq.km")```### Examples of records sharing the same AddressLabel```{python}feature_name ="Addresses_in_fan"gdf = data[feature_name]# Step 1: Find duplicated AddressLabelsduplicates = gdf[gdf.duplicated(subset=["AddressLabel"], keep=False)]# Step 2: Group and inspect (optional)grouped = duplicates.sort_values("AddressLabel").groupby("AddressLabel")# Step 3: Print summaryprint(f"Total duplicated AddressLabel records: {len(duplicates)}")print(f"Number of unique duplicated labels: {duplicates['AddressLabel'].nunique()}")# Step 4: Show some examplesprint("\nExamples of duplicated AddressLabels:")print(duplicates[["AddressLabel", "UnitType", "UnitValue"]].head(10))```### Examples of non-unit addresses with unit parts in AddressLabel```{python}feature_name ="Addresses_in_fan"gdf = data[feature_name].copy()# Step 1: Build a list of distinct, non-null UnitTypespotential_unit_types = gdf["UnitType"].dropna().astype(str).str.strip().unique().tolist()potential_unit_types = [ut for ut in potential_unit_types if ut] # Remove blanks# Step 2: Build regex pattern to match unit suffix at end of address#unit_pattern = r"\b(" + "|".join(map(re.escape, potential_unit_types)) + r")\s+\S+$"unit_pattern =r"\b(?:{})\s+\S+$".format("|".join(map(re.escape, potential_unit_types)))# Step 3: Find rows with UnitType == None and AddressLabel contains a recognizable UnitTypemask = gdf["UnitType"].isnull() & gdf["AddressLabel"].astype(str).str.contains(unit_pattern, regex=True, na=False)# Subset matching rowscandidates = gdf[mask].copy()# Step 4: Reportprint(f"🔍 Found {len(candidates)} rows where AddressLabel ends with a recognizable UnitType + value and UnitType is missing.")print("\n🧪 Sample rows:")print(candidates[["AddressLabel", "UnitType", "UnitValue"]].head())```### Fix these issues found above.```{python}# Step 5: Extract UnitType and UnitValue using regex# Example match: "501 N Arthur Ashe Blvd Unit 2" → ("Unit", "2")unit_extract_pattern =r"\b("+"|".join(map(re.escape, potential_unit_types)) +r")\s+(\S+)$"unit_parts = candidates["AddressLabel"].astype(str).str.extract(unit_extract_pattern)# Assign back to the original gdf using index alignmentgdf.loc[candidates.index, "UnitType"] = unit_parts[0]gdf.loc[candidates.index, "UnitValue"] = unit_parts[1]# Step 6: Update data structuredata[feature_name] = gdfprint(f"\n✅ Extracted UnitType and UnitValue for {len(unit_parts.dropna())} rows.")```### Examples of non-unit addreses sharing the same long/lat```{python}feature_name ="Addresses_in_fan"gdf = data[feature_name]# Step 1: Filter for base addresses (UnitType is None) and valid coordinatesbase_addrs = gdf[ gdf["UnitType"].isnull() & gdf["Longitude"].notnull() & gdf["Latitude"].notnull()].copy()# Step 2: Find duplicate coordinate recordsdup_coords = base_addrs[ base_addrs.duplicated(subset=["Longitude", "Latitude"], keep=False)]# Step 3: Group duplicates for inspectioncoord_groups = dup_coords.groupby(["Latitude", "Longitude"])# Step 4: Summary countsprint(f"📍 Total base address rows with duplicate coordinates: {len(dup_coords)}")print(f"🔁 Unique (lat, lon) pairs with duplicates: {coord_groups.ngroups}")# Step 5: Show sample groupsprint("\n🔍 Example duplicates (UnitType is None, by lat/lon):")for (lat, lon), group in coord_groups:print(f"\nCoordinates: ({lat}, {lon}) — {len(group)} records") sorted_group = group.sort_values("AddressLabel")print(sorted_group[["AddressLabel", "UnitType", "UnitValue", "Latitude", "Longitude"]])break# Remove this if you want to show all groups``````{python}if0:print("\n\n\n🔍 First record from each duplicate (UnitType is None, by lat/lon):")for (lat, lon), group in coord_groups: first = group.sort_values("AddressLabel").iloc[0]print(f"{first['AddressLabel']} | UnitType={first['UnitType']} | UnitValue={first['UnitValue']} | ({lat}, {lon})")print("\n🔍 First record from each duplicate (UnitType is None, by lat/lon):")for (lat, lon), group in coord_groups: first = group.sort_values("AddressLabel").iloc[0] group_size =len(group)print(f"{first['AddressLabel']} | {group_size} records | ({lat}, {lon})")```### Fix: jitter addresses```{python}# Make a copy so we can safely modify coordinatesfrom shapely.geometry import Pointgdf_jittered = gdf.copy()# Define base jitter radius in degrees (very small ~meters)base_radius_deg =0.00001# Loop over each groupfor (lat, lon), group in coord_groups: count =len(group)if count ==1:continue# No need to jitter radius = base_radius_deg * np.sqrt(1.0* count) # Increase radius with count# One random multiplier [0, 1) for each point random_scalars = np.random.rand(count)# Final radius for each point radii = radius * random_scalars angles = np.linspace(0, 2* np.pi, count, endpoint=False) jittered_lats = lat + radii * np.sin(angles) jittered_lons = lon + radii * np.cos(angles)# Assign jittered coordinates back using index gdf_jittered.loc[group.index, "Latitude"] = jittered_lats gdf_jittered.loc[group.index, "Longitude"] = jittered_lonsdata[feature_name].loc[gdf_jittered.index, "Latitude"] = gdf_jittered["Latitude"]data[feature_name].loc[gdf_jittered.index, "Longitude"] = gdf_jittered["Longitude"]# Update all the geometriesdata[feature_name]["geometry"] = data[feature_name].apply(lambda row: Point(row["Longitude"], row["Latitude"]), axis=1)print("✅ Jittered coordinates applied to duplicate lat/lon groups.")```### Example: Unit addresses without non-unit base address```{python}feature_name ="Addresses_in_fan"gdf = data[feature_name]# Step 1: Get unit addresses (where UnitType is not null)unit_addrs = gdf[gdf["UnitType"].notnull()].copy()# Step 2: Extract base address (remove unit from AddressLabel)# Assumes unit info appears after a comma, e.g., "123 Main St, Apt 2B"unit_addrs["BaseAddressLabel"] = ( unit_addrs["BuildingNumber"].fillna("") +" "+ unit_addrs["StreetDirection"].fillna("") +" "+ unit_addrs["StreetName"].fillna("") +" "+ unit_addrs["StreetType"].fillna("")).str.replace(r"\s+", " ", regex=True).str.strip()# Step 3: Get base addresses in the dataset (with UnitType == None)base_addrs = gdf[gdf["UnitType"].isnull()].copy()base_addrs["BaseAddressLabel"] = ( base_addrs["BuildingNumber"].fillna("") +" "+ base_addrs["StreetDirection"].fillna("") +" "+ base_addrs["StreetName"].fillna("") +" "+ base_addrs["StreetType"].fillna("")).str.strip().replace(r"\s+", " ", regex=True)base_addrs["AddressLabel"] = base_addrs["AddressLabel"].fillna("").astype(str).str.strip()# Step 4: Get the set of base labels from base_addrsbase_set =set(base_addrs["BaseAddressLabel"])# Step 5: Find unit addresses whose base is missingmissing = unit_addrs[~unit_addrs["BaseAddressLabel"].isin(base_set)]# Step 6: Print summaryprint(f"Total unit addresses: {len(unit_addrs)}")print(f"Unit addresses with missing base address: {len(missing)}")ifnot missing.empty:print("\nExamples of missing base addresses:")print(missing[["AddressLabel", "BaseAddressLabel", "UnitType", "UnitValue"]].head())```### List of non-unit addresses where baseaddresslabel <> addressLabel```{python}# Count total base recordstotal_base =len(base_addrs)# Compare BaseAddressLabel to AddressLabelmatch = base_addrs[base_addrs["BaseAddressLabel"] == base_addrs["AddressLabel"]]mismatch = base_addrs[base_addrs["BaseAddressLabel"] != base_addrs["AddressLabel"]]# Print countsprint(f"Total base addresses (UnitType is null): {total_base}")print(f"Matching AddressLabel == BaseAddressLabel: {len(match)}")print(f"Non-matching AddressLabel != BaseAddressLabel: {len(mismatch)}")# Show examples of mismatchesifnot mismatch.empty:print("\nExamples of mismatched base addresses:")print(mismatch[["AddressLabel", "BaseAddressLabel"]].head())```### Addresses with the top 20 unit counts.```{python}feature_name ="Addresses_in_fan"gdf = data[feature_name].copy()# Step 1: Create BaseAddressLabel for all addressesgdf["BaseAddressLabel"] = ( gdf["BuildingNumber"].fillna("") +" "+ gdf["StreetDirection"].fillna("") +" "+ gdf["StreetName"].fillna("") +" "+ gdf["StreetType"].fillna("")).str.strip().replace(r'\s+', ' ', regex=True)# Step 2: Count number of unit addresses per BaseAddressLabelunit_counts = ( gdf[gdf["UnitType"].notnull()] .groupby("BaseAddressLabel") .size() .rename("UnitCountForBase"))# Step 3: Map unit counts to base addressesgdf["UnitCount"] = gdf["BaseAddressLabel"].map(unit_counts).fillna(1).astype(int)# Save back to data objectdata[feature_name] = gdf``````{python}feature_name ="Addresses_in_fan"gdf = data[feature_name]# Sort by UnitCount descending and drop duplicates to avoid listing each unit individuallytop_bases = ( gdf[gdf["UnitType"].isnull()] # Only base addresses .sort_values("UnitCount", ascending=False) .head(20))# Display relevant infoprint(top_bases[["AddressLabel", "BaseAddressLabel", "UnitCount","UnitType"]])```## FDA - Static MapBelow is a static map of the Fan District Association, using dataobtained from the City of Richmond GeoHub. The image shows the boundaryof the Fan District Association in black. Addresses are identifed by red dots.Parcels are identified as gray regions.[Click here](./fan_map.png) to download a PNG of this image.::: {.column-page-inset-right}```{python}# Plottingfig, ax = plt.subplots(figsize=(10, 10))data[selector_key].boundary.plot(ax=ax, color="black", linewidth=3, label="The Fan boundary")data["Parcels_in_fan"].plot(ax=ax, color="lightgray", edgecolor="gray", alpha=0.7, label="Parcels")data["Addresses_in_fan"].plot(ax=ax, color="red", markersize=5, label="Addresses")#data["Neighborhoods_in_fan"].plot(ax=ax, facecolor='none', edgecolor="lightblue", markersize=5, label="Neighborhoods")legend_elements = [ Line2D([0], [0], color="black", lw=3, label="The Fan boundary"), Patch(facecolor="lightgray", edgecolor="gray", label="Parcels"), Line2D([0], [0], marker='o', color='w', label="Addresses", markerfacecolor='red', markersize=6),# Line2D([0], [0], color="lightblue", lw=1.5, label="Neighborhoods boundary")]ax.legend(handles=legend_elements)ax.set_title("Addresses and Parcels in The Fan")ax.axis("off")plt.tight_layout()ax.set_title("Addresses and Parcels in The Fan")ax.axis("off")plt.tight_layout()plt.savefig("../docs/fan_map.png", dpi=300, bbox_inches="tight")plt.show()```:::## Interactive Map: Addresses and Parcels in The FanBelow is an interactive map of the Fan District Association. You canzoom and scroll.```{python}import foliumfrom shapely.geometry import mappingimport geopandas as gpdif0:# Create a base map centered on The Fan projected = data[selector_key].to_crs(epsg=2283)# Compute centroid in projected CRS (accurate) centroid = projected.geometry.centroid.iloc[0]# Convert back to geographic CRS (WGS84) for Folium centroid_latlon = gpd.GeoSeries([centroid], crs=2283).to_crs(epsg=4326).geometry.iloc[0] fan_center = centroid_latlon.coords[0][::-1] # (lat, lon)# Create map m = folium.Map(location=fan_center, zoom_start=15, tiles="cartodbpositron")# Project and compute bounding box in WGS84fan_shape = data[selector_key].to_crs(epsg=4326) # Folium uses WGS84# Compute bounds: [[south, west], [north, east]]minx, miny, maxx, maxy = fan_shape.total_boundsbounds = [[miny, minx], [maxy, maxx]]# Center for initial rendering (optional fallback)center = [(miny + maxy) /2, (minx + maxx) /2]# Create map and set boundsm = folium.Map(location=center, zoom_start=15, tiles="cartodbpositron")m.fit_bounds(bounds)``````{python}# Add neighborhood boundaryx = folium.GeoJson( data[selector_key].geometry, name="The Fan Boundary", style_function=lambda x: {"color": "black","weight": 3,"fillOpacity": 0, }).add_to(m)``````{python}# Add parcels (lighter gray polygons)x = folium.GeoJson( data["Parcels_in_fan"].geometry, name="Parcels", style_function=lambda x: {"color": "#999999","weight": 0.5,"fillOpacity": 0.4, },).add_to(m)# Add addresses (red points)``````{python}from shapely.geometry import Pointdef get_color_for_unit_count(unit_count: int, max_unit_count: int=10) ->str:if unit_count ==1:return"#FF0000"# red# For unit_count > 1: magenta to white gradient scale =min(unit_count, max_unit_count) / max_unit_count green_value =int(255* scale) # from 0 (magenta) to 255 (white)returnf"#FF{green_value:02X}FF"for _, row in data["Addresses_in_fan"].iterrows(): pt = row.geometry# Skip labels for addresses where UnitType is not None. unit_type = row["UnitType"]if unit_type isnotNone:continue# Only plot if geometry is a valid Pointifnotisinstance(pt, Point):print(f"Skipping non-Point geometry at index {_}: {type(pt)}")continueif (row.get("BuildingNumber")=="1465") and (row.get("StreetName")=="Floyd"):#logger.debug( row )pass unit_count = row.get("UnitCount",1) label = row.get("AddressLabel", "")if unit_count>1: label = label +f" ({unit_count})" tooltip = folium.Tooltip(label) if label.strip() elseNone color = get_color_for_unit_count(unit_count) folium.CircleMarker( location=[pt.y, pt.x], radius=2, color=color, fill=True, fill_opacity=0.8, tooltip=tooltip, # Explicit safe wrapper ).add_to(m)# Add layer control and display mapx = folium.LayerControl().add_to(m)```::: {.column-page-inset-right}```{python}#| fig-height: 10#| fig-width: 12m```## Interactive Map: Heatmap of Addresses and Unitcounts```{python}# Project and compute bounding box in WGS84fan_shape = data[selector_key].to_crs(epsg=4326) # Folium uses WGS84# Compute bounds: [[south, west], [north, east]]minx, miny, maxx, maxy = fan_shape.total_boundsbounds = [[miny, minx], [maxy, maxx]]# Center for initial rendering (optional fallback)center = [(miny + maxy) /2, (minx + maxx) /2]# Create map and set boundsm = folium.Map(location=center, zoom_start=15, tiles="cartodbpositron")m.fit_bounds(bounds)# Add neighborhood boundaryx = folium.GeoJson( data[selector_key].geometry, name="The Fan Boundary", style_function=lambda x: {"color": "black","weight": 3,"fillOpacity": 0, }).add_to(m)# Add parcels (lighter gray polygons)x = folium.GeoJson( data["Parcels_in_fan"].geometry, name="Parcels", style_function=lambda x: {"color": "#999999","weight": 0.5,"fillOpacity": 0.4, },).add_to(m)# Add intensities for heatmap layerfrom folium.plugins import HeatMap# Step 1: Project to WGS84 (if not already)addresses = data["Addresses_in_fan"].to_crs(epsg=4326)# Step 2: Build heatmap data: [[lat, lon, intensity], ...]# Clip UnitCount at 150 to suppress outlier skewclip_value =90addresses["ClippedUnitCount"] = addresses["UnitCount"].clip(upper=clip_value)# Normalizeheat_data = [ [geom.y, geom.x, row.ClippedUnitCount / clip_value]for geom, row inzip(addresses.geometry, addresses.itertuples())if geom isnotNoneand row.ClippedUnitCount >0and row.UnitType isNone]# Step 3: Add HeatMap layerHeatMap( heat_data, radius=12, blur=10, max_zoom=15, name="High-Density Housing").add_to(m)m``````{python}#m.save("../docs/fan_map.html")#m.save("./fan_map.html")#m```:::<!-- see: https://quarto.org/docs/authoring/article-layout.html --><!--::: {.column-page-inset-right}<iframe src="./fan_map.html" width="100%" height="750px" style="border:none;" data-external="1" ></iframe>:::-->## Interactive Map: Parcels and UnitCounts```{python}# Step 1: Filter addresses with UnitType == Noneaddresses = data["Addresses_in_fan"]addresses = addresses[addresses["UnitType"].isna()].copy()# Ensure both are in the same CRSaddresses = addresses.to_crs(data["Parcels_in_fan"].crs)# Step 2: Perform spatial join - match addresses within parcelsjoined = gpd.sjoin(addresses, data["Parcels_in_fan"], how="inner", predicate="within")# Step 3: Group by Parcel ID (or index) and sum UnitCount# Choose a unique key to group by — use index_right if no parcel ID is availableparcel_unit_counts = ( joined.groupby("index_right")["UnitCount"] .sum() .rename("SummedUnitCount"))# Step 4: Join the result back to Parcels GeoDataFrameparcels_with_units = data["Parcels_in_fan"].copy()parcels_with_units["SummedUnitCount"] = parcel_unit_counts# Optional: fill NaN with 0 for parcels without addressesparcels_with_units["SummedUnitCount"] = parcels_with_units["SummedUnitCount"].fillna(0).astype(int)``````{python}fan_shape = data[selector_key].to_crs(epsg=4326) # Folium uses WGS84# Compute bounds: [[south, west], [north, east]]minx, miny, maxx, maxy = fan_shape.total_boundsbounds = [[miny, minx], [maxy, maxx]]# Center for initial rendering (optional fallback)center = [(miny + maxy) /2, (minx + maxx) /2]# Create map and set boundsm = folium.Map(location=center, zoom_start=15, tiles="cartodbpositron")m.fit_bounds(bounds)# Add neighborhood boundaryx = folium.GeoJson( data[selector_key].geometry, name="The Fan Boundary", style_function=lambda x: {"color": "black","weight": 3,"fillOpacity": 0, }).add_to(m)# Add parcels (lighter gray polygons)x = folium.GeoJson( data["Parcels_in_fan"].geometry, name="Parcels", style_function=lambda x: {"color": "#999999","weight": 0.5,"fillOpacity": 0.4, },).add_to(m)if0:# Step 1: Clip values (saturate above 100) clip_max =100 parcels_with_units["ClippedUnits"] = parcels_with_units["SummedUnitCount"].clip(upper=clip_max)# Step 2: Create colormap from 0 to clip_max clipped_colormap = cm.linear.YlOrRd_09.scale(0, clip_max) clipped_colormap.caption =f"Units per Parcel (capped at {clip_max})" clipped_colormap.add_to(m)# Step 3: Add to map folium.GeoJson( parcels_with_units, name=f"Units per Parcel (max {clip_max})", style_function=lambda feature: {"fillColor": clipped_colormap(feature["properties"]["ClippedUnits"])if feature["properties"]["ClippedUnits"] else"#ffffff","color": "#666666","weight": 0.5,"fillOpacity": 0.7, },# tooltip=folium.GeoJsonTooltip(# fields=["SummedUnitCount", "ClippedUnits"],# aliases=["Total Units", f"Clipped to {clip_max}"],# localize=True# ) ).add_to(m)if0:import branca.colormap as cm colormap = cm.linear.YlOrRd_09.scale(0, parcels_with_units["SummedUnitCount"].max()) colormap.caption ="Units per Parcel (filtered by UnitType=None)" colormap.add_to(m) x = folium.GeoJson( parcels_with_units, name="Filtered Unit Density", style_function=lambda feature: {"fillColor": colormap(feature["properties"]["SummedUnitCount"]) if feature["properties"]["SummedUnitCount"] else"#ffffff","color": "#333333","weight": 0.5,"fillOpacity": 0.7, }, tooltip=folium.GeoJsonTooltip( fields=["ParcelID", "SummedUnitCount"], aliases=["Parcel", "Total Units"], ) ).add_to(m)if1:def classify_unit_bin(unit_count):if unit_count ==0:return'0'elif unit_count ==1:return'1'elif unit_count ==2:return'2'elif unit_count <=4:return'3-4'elif unit_count <=10:return'5-10'elif unit_count <=20:return'11-20'elif unit_count <=35:return'21-35'else:return'36+'# Apply to DataFrame parcels_with_units["UnitBin"] = parcels_with_units["SummedUnitCount"].apply(classify_unit_bin)# format dollar values parcels_with_units["FormattedValue"] = parcels_with_units["TotalValue"].apply( lambda x: f"${int(x/1000)}k"if x else"$0k" )# Define bin colors bin_colors = {'0': '#ffffff', # white'1': '#ffffcc', # light yellow'2': '#fd8d3c', # orange'3-4': '#fc4e2a', # darker orange'5-10': '#e31a1c', # red'11-20': '#bd0026', # dark red'21-35': '#800026', # deeper red'36+': '#4d0018', # very dark red } folium.GeoJson( parcels_with_units, name="Binned Unit Density", style_function=lambda feature: {"fillColor": bin_colors.get(feature["properties"]["UnitBin"], "#cccccc"),"color": "#666666","weight": 0.5,"fillOpacity": 0.7, }, tooltip=folium.GeoJsonTooltip( fields=["ParcelID","SummedUnitCount", "UnitBin"], aliases=["ParcelID","Total Units", "Bin"], localize=False ) ).add_to(m)m```## Available columns```{python}#| output: asisprint(":::: {.columns} ")width =round(floor(100.0/len( [selector] + features )))for feature in [selector] + features:if feature in features: feature = feature +"_in_fan"print(f"""::: {{.column width={width}%}}### {feature}| Property ||----------|""") columns = data[feature].columns.tolist()for col in columns:print(f"| {col} |")print (f"""::: """)print("\n::::")```
